import pandas as pd

import numpy as np

import random

import matplotlib.pyplot as plt

import os

import time


# 定义模拟退火算法类

class SimulatedAnnealing:

    def __init__(self, initial_params, param_ranges, max_iterations, initial_temp, cooling_rate, data):

        self.params = initial_params  # 初始参数

        self.param_ranges = param_ranges  # 参数范围

        self.max_iterations = max_iterations  # 最大迭代次数

        self.temp = initial_temp  # 初始温度

        self.cooling_rate = cooling_rate  # 降温速率

        self.data = data  # 数据

        self.best_solution = initial_params.copy()  # 最优解

        self.best_error, self.best_mae, self.best_rmse = self.objective_function(**initial_params)  # 最优误差

        self.error_history = []  # 记录每次迭代的误差


    # 目标函数

    def objective_function(self, k, a, b, c1, c2):

        P_actual = self.data['P']

        f = self.data['f']

        B = self.data['B']

        t = self.data['t']

        P_predicted = k * (f ** a) * (B ** b) * (1 + c1 * t + c2 * (t ** 2))


        # 相对误差

        relative_error = np.mean(((P_actual - P_predicted) / P_actual) ** 2)

        # 平均绝对误差 (MAE)

        mae = np.mean(np.abs(P_actual - P_predicted))

        # 均方根误差 (RMSE)

        rmse = np.sqrt(np.mean((P_actual - P_predicted) ** 2))


        return relative_error, mae, rmse


    # 执行一次模拟退火迭代

    def run(self):

        for iteration in range(self.max_iterations):

            # 生成新的解

            new_params = self.generate_new_solution(self.params)

            new_error, new_mae, new_rmse = self.objective_function(**new_params)


            # 比较新解与当前解

            if new_error < self.best_error:

                self.params = new_params

                self.best_error = new_error

                self.best_mae = new_mae

                self.best_rmse = new_rmse

                self.best_solution = new_params.copy()

            else:

                # 按一定概率接受较差解

                if self.accept_solution(new_error):

                    self.params = new_params


            # 降低温度

            self.temp *= self.cooling_rate

            self.error_history.append(self.best_error)


        return self.best_solution, self.best_error, self.best_mae, self.best_rmse


    # 生成新解

    def generate_new_solution(self, current_params):

        new_params = current_params.copy()

        for param_name, (min_val, max_val) in self.param_ranges.items():

            perturbation = random.uniform(-0.1, 0.1)

            new_params[param_name] += perturbation

            new_params[param_name] = max(min_val, min(max_val, new_params[param_name]))


        return new_params


    # 接受新解的概率

    def accept_solution(self, new_error):

        delta_error = new_error - self.best_error

        acceptance_probability = np.exp(-delta_error / self.temp)

        return random.random() < acceptance_probability


# 记录运行开始时间

start_time = time.time()


# 读取CSV文件路径

file_path = r'D:\2024华为杯数学建模\2024C题-数据驱动与电磁损耗\C题\附件一数据拆解\Q2\Q2磁芯材料1正弦波形.csv'

if not os.path.exists(file_path):

    raise FileNotFoundError(f"文件不存在: {file_path}")

try:

    data = pd.read_csv(file_path, skiprows=1, encoding='gbk')

except UnicodeDecodeError:

    data = pd.read_csv(file_path, skiprows=1, encoding='ISO-8859-1')


data.columns = ['t', 'f', 'P', 'Unknown', 'B']  # 重命名列


# 定义参数范围

param_ranges = {

    'k': (0.05, 0.15),  # k 的范围

    'a': (1.5, 1.8),  # a 的范围

    'b': (2.5, 2.7),  # b 的范围

    'c1': (-1, 1),  # c1 的范围

    'c2': (-1, 1)  # c2 的范围

}


# 初始参数

initial_params = {

    'k': random.uniform(param_ranges['k'][0], param_ranges['k'][1]),  # k 的初始值

    'a': random.uniform(param_ranges['a'][0], param_ranges['a'][1]),  # a 的初始值

    'b': random.uniform(param_ranges['b'][0], param_ranges['b'][1]),  # b 的初始值

    'c1': random.uniform(param_ranges['c1'][0], param_ranges['c1'][1]),  # c1 的初始值

    'c2': random.uniform(param_ranges['c2'][0], param_ranges['c2'][1])  # c2 的初始值

}


# 保存每次迭代结果

results = []


# 创建输出目录

output_dir = r'D:\2024华为杯数学建模\2024C题-数据驱动与电磁损耗\C_code\Q2_code\Q2可视化结果\Q2修正图示\二次函数\TU'

os.makedirs(output_dir, exist_ok=True)


# 运行50次模拟退火算法

for run in range(50):

    # 实例化模拟退火算法

    sa = SimulatedAnnealing(

        initial_params=initial_params,

        param_ranges=param_ranges,

        max_iterations=1000,  # 设置迭代次数

        initial_temp=1000,  # 初始温度

        cooling_rate=0.95,  # 降温速率

        data=data

    )


    # 运行模拟退火算法

    best_solution, best_error, best_mae, best_rmse = sa.run()


    # 将每次结果保存到列表

    results.append([run + 1, best_solution['k'], best_solution['a'], best_solution['b'], best_solution['c1'], best_solution['c2'], best_error, best_mae, best_rmse])


    print(f"第{run + 1}次运行结果: k = {best_solution['k']}, a = {best_solution['a']}, b = {best_solution['b']}, c1 = {best_solution['c1']}, c2 = {best_solution['c2']}, 最小误差 = {best_error}, MAE = {best_mae}, RMSE = {best_rmse}")


    # 可视化实际值与预测值

    f = data['f']

    B = data['B']

    t = data['t']

    P_actual = data['P']

    P_predicted = best_solution['k'] * (f ** best_solution['a']) * (B ** best_solution['b']) * (1 + best_solution['c1'] * t + best_solution['c2'] * (t ** 2))


    plt.figure(figsize=(10, 6))

    plt.scatter(P_actual, P_predicted, color='blue', label='preP')

    plt.plot([P_actual.min(), P_actual.max()], [P_actual.min(), P_actual.max()], 'r--', label='y=x')

    plt.xlabel('trueP')

    plt.ylabel('preP')

    plt.title(f'prePVStrueP - run_num {run + 1}')

    plt.legend()

    plt.grid()

    plt_path = os.path.join(output_dir, f'真实P_vs_预测P_{run + 1}.png')

    plt.savefig(plt_path)

    plt.close()


    # 绘制迭代曲线

    plt.figure(figsize=(10, 6))

    plt.plot(range(1, len(sa.error_history) + 1), sa.error_history, label='best error')

    plt.xlabel('iteration')

    plt.ylabel('error')

    plt.title(f'error evolution {run + 1}')

    plt.legend()

    plt.grid()

    convergence_path = os.path.join(output_dir, f'迭代曲线_运行_{run + 1}.png')

    plt.savefig(convergence_path)

    plt.close()


# 记录运行结束时间

end_time = time.time()


# 计算并输出运行时间

execution_time = end_time - start_time

print(f"运行50次的总时间: {execution_time:.2f} 秒")


# 将所有运行结果保存到 CSV 文件

output_csv_path = os.path.join(output_dir, '50次拟合结果.csv')

result_df = pd.DataFrame(results, columns=['run_num', 'k', 'a', 'b', 'c1', 'c2', 'relative_error', 'mae', 'rmse'])

result_df.to_csv(output_csv_path, index=False, encoding='utf-8-sig')

print(f"所有50次拟合结果已保存到: {output_csv_path}")